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We present an empirical strategy to determine the Hamiltonian dynamics of a two-qubit system us- 
ing only initialization and measurement in a single fixed basis. Signal parameters are estimated from 
measurement data using Bayesian methods from which the underlying Hamiltonian is reconstructed, 
up to three unobservable phase factors. We extend the method to achieve full control Hamiltonian 
tomography for controllable systems via a multi-step approach. The technique is demonstrated and 
evaluated by analyzing data from simulated experiments including projection noise. 
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I. INTRODUCTION 

Using quantum phenomena to perform new modes 
of computation is a daunting challenge [T]. Signifi- 
cant achievements in the theory of quantum computa- 
tion include the development of error correction, fault- 
tolerance [2], and scalability of quantum circuits [3]. 
However, in order to build large scale quantum proces- 
sors, many individual quantum systems must be manip- 
ulated with extraordinary precision and accuracy. A pre- 
requisite for this level of quantum control is precise char- 
acterization of the underlying dynamics and its response 
to control fields, so-called Hamiltonian Engineering (|H|S] 
and references therein). This is especially crucial for 
manufactured devices such as solid state quantum bits 
(qubits), e.g. quantum dots (Fig. [ij or superconducting 
quantum interference devices (SQUIDs). Any manufac- 
turing process will introduce variations so it is important 
to empirically identify the control relationship for each 
component. In a large-scale quantum computer, it is de- 
sirable to be able to achieve this using in situ resources, 
i.e., initialization, control actuators and measurement ca- 
pabilities already present for performing computation. 

The canonical method for assessing quantum dynamics 
is quantum process tomography (QPT) |51[71|5]. This in- 
volves initialization of a quantum system in a (complete) 
set of states, allowing it to evolve under the dynamics un- 
der consideration, and then performing an information- 
ally complete measurement on the output state for each 
input. From this set of input-output data, the superop- 
erator, or completely positive (CP) map, governing the 
quantum evolution of the system can be reconstructed. 
This may then be repeated for different evolution times 
to obtain an estimate of the Lindblad operators (gener- 
ators of the dynamics) [9]. For control purposes, QPT 
would be performed for a variety of actuator settings to 
build up a map of the control space. 
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FIG. 1: Manufactured Qubit System. A pair of horizontally 
aligned double quantum dots (center) can act as a two-qubits. 
A qubit can be defined in each double quantum dot by two 
different charging states, e.g. a single excess electron located 
on the left or right dot of each pair. Electrodes (top and 
bottom) control the potentials and electron tunneling rates. 
Single electron transistors (left and right) measure the loca- 
tions of the excess electrons which defines the measurement 
basis, or logical states of the qubit. Due to finite manufac- 
turing precision, the placement of the control and measure- 
ment structures may not be exactly as calculated, hence the 
Hamiltonian dependence upon control signals will have to be 
determined empirically. Image courtesy of Hitachi Cambridge 
Laboratories, Hitachi Europe Ltd. 



A potential disadvantage of QPT is the need for ah 
initio initialization and measurement outside of the com- 
putational basis, a capability which may not exist in the 
absence of characterization in the first place. It is usu- 
ally argued that initialization and measurement in an 
arbitrary basis can be achieved by unitary rotation of a 
fixed basis, however this pre-supposes that the system 
response to control fields has already been characterized, 
a vicious circle. Previous work has addressed this issue 
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for the case of a single qubit subject to multiple control 
Hamiltonians, decoherence, and imperfect subspace con- 
finement [ini EH m m [11 [H]. Here, we extend the 
basic idea of Hamiltonian characterization to two cou- 
pled qubits with an unknown generic internal Hamilto- 
nian and control Hamiltonian response. 

This paper is organized as follows. In Section |IT] we 
discuss the basic principles of Hamiltonian tomography 
for a two-qubit system with a fixed but unknown Hamil- 
tonian, assuming only the ability to measure the system 
at specific times in a fixed measurement basis, but no 
control or a priori knowledge of the system. We also 
deliberately exclude the ability to perform local opera- 
tors on either qubit, or the ability to initialize the sys- 
tem in states other than the measurement basis states. 
Our approach differs in this regard from related work on 
two-qubit Hamiltonian identification using concurrence 
spectroscopy [HI [T^ or optimal experiment design [l8] . 
These approaches may be preferable for certain types of 
systems but have some limitations as they presume the 
single qubit dynamics can be fully characterized inde- 
pendently of the inter-qubit coupling, which is required 
to prepare the two-qubit system in superposition states 
by applying local rotations. Using concurrence also limits 
us to reconstructing the non-local part of the two-qubit 
Hamiltonian. Thus, this approach may be well-suited for 
some systems e.g., with weak-coupling and non-local in- 
teraction Hamiltonians, but may be problematic for other 
systems. The approach in this paper should be seen as 
complementary to these works. 

In Section IIIII we discuss how to extract the relevant 
system parameters from the noisy measurement data, 
accurately and robustly. The difficulty of this task is 
greatly magnified compared to the single qubit case due 
to the number of parameters to be determined, as well as 
the increased signal complexity. A nai've approach using 
straightforward least-squares error minimization failed 
completely when applied to noisy data from simulated ex- 
periments. The power spectrum of the signal (which was 
sufficient for the single qubit case) is still useful, but no 
longer an optimal frequency estimator in the presence of 
multiple frequencies, and obtaining accurate estimates of 
the amplitudes of different frequency components is very 
difficult. For these reasons Bayesian analysis is employed 
to determine the signal parameters, which is shown to 
result in significant improvements in the accuracy and 
robustness of the procedure. 

In Section IIVI we show how to reconstruct the total 
Hamiltonian, or more precisely, its matrix representation 
with respect to the fixed measurement basis, from the 
estimated parameters. Unlike the single qubit case, cal- 
culating the 16 matrix elements of the two-qubit Hamil- 
tonian from the 214 parameters estimated from the 16 
measured signals is non-trivial, and requires several op- 
timization steps, from identifying the most likely level 
structure from the set of transitions frequencies, to de- 
termining the magnitudes and phases of the Hamiltonian 
matrix elements that provide the best fit with the esti- 



mated parameters. The analysis also shows that the fixed 
Hamiltonian can be determined only up to a global phase 
and sign, as well as three phases, which define U(l) trans- 
formations of the measurement basis states. If there are 
no other measurements or control available then these 
U(l) transformations of the basis states have no observ- 
able effect. Modulo these unobservable parameters, we 
demonstrate that we can reconstruct the overall Hamil- 
tonian with very good accuracy from noisy data. 

In Section |V] we consider the more general case of 
control Hamiltonian tomography. In particular, we are 
interested in characterizing Hamiltonians H = ^^(f) 
that depend on a number of external parameters f — 
(/i, . . . , /m) that can be varied experimentally, such as 
voltages applied to certain gate electrodes that allow us 
to vary confinement potentials, tunneling rates etc. By 
varying these parameters over time we can engineer com- 
plicated effective Hamiltonians and efficiently achieve a 
wide range of control tasks from quantum state prepara- 
tion to gate implementation |19j using powerful optimal 
control techniques |20l [21]. However, effective control 
requires knowledge of the dependence of the Hamilto- 
nian on these parameters H(f). When applying different 
Hamiltonians, the previously unobservable phase factors 
now have practical effects and are critical for full control 
Hamiltonian tomography. We show how to determine 
these phases, relative to a reference Hamiltonian, using 
a simple two-step experiment, and how to use this infor- 
mation to achieve full control-Hamiltonian tomography. 

Finally, in Section [VT| we discuss applications of the re- 
sults, as well as future improvements and generalizations 
to our method. 



II. FIXED HAMILTONIAN TOMOGRAPHY 

Throughout this paper we assume that we are given 
a two-qubit system with an unknown Hamiltonian, and 
a measurement apparatus that enables us to perform a 
fixed projective measurement on each qubit, including 
the ability to perform effectively simultaneous measure- 
ments on both qubits [57j. We denote the measurement 
basis states of the resulting four-outcome measurement 
by |1) = |00), |2) = |01), |3>^ = |10) and |4) = |11). We 
then perform the following simple experiment: 

1. Initialize the system in one of the four measure- 
ment basis states \k) by performing simultaneous 
measurements on both qubits. 

2. Let the system evolve for some time t. 

3. Perform simultaneous measurements on both 
qubits, projecting the system back into one of the 
four measurement basis states. 

By repeating this experiment many times for a fixed evo- 
lution time t, we can estimate the probabilities pteit) = 
I (^l^'fe(t)) P, where [^'^(i)) is the time-evolved state and 
l^fc(O)) = \k). By further repeating the experiment for 
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different times t„ for n = 0, . . . , N — 1 we can strobo- 
scopically capture the evolution of the probabiHties Pke{t) 
for k,i — 1,2,3,4, yielding 16 noisy signals as shown in 
Fig. [2] [28 . What information about the Hamiltonian can 
we extract from this data, and what is the most effective 
way to extract this information? 

Assume the evolution of the system is governed by a 
fixed Hamiltonian according to the Schrodinger equation 



th-\nt)) = H\nt)). 



(1) 



Expanding the Hamiltonian H with respect to its or- 
thonormal eigenbasis : ~ 1, . . . , 4}, 



(2) 



where Xi, are the (real) eigenvalues, and setting 



rkuB 



we obtain 



{e\H\k) = ^ A,(£|C.)(e.|fc) = ^ A.r.^rfc.e^^*--'^'''). 



i/=i 



Further defining Ske;u = ruuTi^, 5m-u = 



(3) 

hv " 4>i>i> and 



^ki-nu = -Ski;ii + 5kt,v we obtain 



{t\H\k) ^Y.^-'ki-ue 



JSi, 



where the phase terms satisfy 5(k:i = —^kt.\ and 



(4) 



23;1 


= ^13;1 


- <5l2;l, 


(5a) 


24;1 


= ^14;1 


- <5l2;l, 


(5b) 


34; 1 


= ^14;1 


- <5l3;l. 


(5c) 



If the system is initialized in one of the measurement 
basis states |^'fc(0)) = |/c), its time-evolved state \^k{i)) 
under the action of H is given by 



|vl/,(t))=^e-^^*|U(e.|fc) 



(6) 



and since (^,^1$;^) = ^vji, its projection onto the measure- 
ment basis state \C) at time t is 



(7) 



Hence, the probability Pke{t) — |(^|^'fc(t))P of the out- 
come \£) for a projective measurement of |^'fe(t)) is 



E 



Ski-l'S 



E ^kt-fiB 



E + ^ E ^kl:^lSke■^'COs{u!^,yt - Afc£;^,,), 



where w^i^ 



A. 



Ap, and using cos (a — b) 



cos(a) cos(6) + sin(a) sin(6), 
Pki{t) = Cki + 2 E afe^;^!^ cos(t^^i.t) -f- bkt,t^„ sm{uj^,yt) 

(8) 



where the coefficients are 



o/cf:^!/ = Ske;uSke;fj.cos{Ake-i_t,^) (9a) 
= Sfc£;^Sfef;i.sin(Afe£;p^) (9b) 

Cki = J2^4e;u- (9c) 

Eq. (|9| shows that the observed dynamics are com- 
pletely determined by the transition frequencies ojfj,^, the 
phase differences Ake-p,i, and the (real) coefficients akt.^m, 
bke-fif and Cki, from which we can reconstruct the Hamil- 
tonian H defined by 



{e\H\k) = Y,\ 



(10) 



u=l 



where Xi, — ujiu — \ {loi2 +^^13 +^14), which is related to 
the actual Hamiltonian H by 

H = D^HD + const. I, (11) 

where D = diag(l, e^^i" , e^^i^ ^ e*"^!-*). The last term is sim- 
ply a global energy shift which has no observable conse- 
quences in general. The diagonal operator D represents 
the U(l) degree of freedom for redefining each measure- 
ment basis state. With only a single constant Hamilto- 
nian, and preparation and measurement in a single fixed 
basis only, we cannot completely determine the Hamilto- 
nian. 



III. PARAMETER ESTIMATION 

The first task is to analyze the measurement traces 
Eq. ([S]) and extract signal parameters Eq. ^ and the 
frequencies w^i/. For convenience we label the transition 
frequencies of the system Wm for m = 1, . . . , 6, assuming 
ujm+i > > 0, and define the vectors uj = (wi, . . . , loq), 
a/cf = (afcf;m) and hke = (6fc£;m) for fc,^ = 0, 1,2,3. The 
first step towards identifying the Hamiltonian H is to 
extract the six transition frequencies a; and 13 linear co- 
efficients afe£, hki-^m, and Cki for each of the 16 signals. 
Although there are 6 -|- 13 x 16 = 214 parameters, the 
problem would be relatively simple if pki {t) was known 
with infinite precision for a set of sample times t„. In 
practice, the accuracy of Pke{tn) is limited by noise, in 
our case projection noise due to the finite number of rep- 
etitions Ne, which renders the problem one of parameter 
estimation for a harmonic signal with multiple frequen- 
cies and phases. Problems of this type are common in en- 
gineering from acoustics to image processing, and many 
techniques have been developed, but our parameter esti- 
mation problem is non-trivial due to the large number of 
parameters involved. 
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FIG. 2: Simulated measurements of System 1 with 2+1 data points per trace sampled at At — 0.1, signal length T = 102.4 
(arbitrary units), number of experiment repetitions per data point A^e = 250. Each graph is the probability pki{t) at time t of 
detecting the system in state \l) {i = 1, 2, 3, 4 left to right) if initialized in state |fe) (fe = 1, 2, 3, 4 top to bottom). 



According to Eq. ^ the traces Pkeit) should be lin- 
ear combinations of the 13 basis functions g2m-i(i) = 
cos{ujmt), g2m(t) = sin(wmi) for TO = 1,...,6, and 
gisit) ^ 1, i.e., 



ake,mg2m-l{t) + bk£,mg2m{t) + Cki (12) 



and our objective is to find parameters ujm, cike-m, bki-m 
and Cki that maximize the likelihood of the measured 
data. Setting d^^ = (dke-i, ■ ■ ■ ,dki-N), where dke-n de- 
notes the approximate value of Pkt.{tn) derived from the 
measurement data, one way to proceed is to try to fit the 
parameters to minimize the squared L^-norm of the error 



k,e. 



|efc£ 



El 



\Pkl 



Ike 



|2 
|2j 



(13) 



where pki = (pke-i, ■ ■ • ,Pkt,N) with pkt-n = Pkt{tn) and 
ll^lli — Ylin=i^n as usual. However, for problems with 
a large number of noisy data points and a large number 
of parameters, as in our case, finding a solution close to 
the (unknown) global minimum of the error using brute- 
force optimization over all system parameters at once is 
difficult at best. We tested this strategy and in most 
cases achieved only poor results. 

Instead of minimizing the global error, we can alterna- 



tively try to maximize the related likelihood function 



L{aiki,hkt,Cki,ijJ,(j) ^ W crj.^^exp 



kj=l 



iPke 



(14) 

Note that we have implicitly assumed here that the 
signals Pke{t) are independent and subject to Gaussian 
white noise with variance cr^^, assumptions that are not 
strictly valid in our case. Hermitian symmetry of the 
Hamiltonian requires a^^ = a^j. and hke = ~^ik but we 
will enforce this symmetry later by averaging the esti- 
mated coefficients 



{hki — hik)- 



(15a) 
(15b) 



The Gaussian noise model is not strictly valid either; 
if the measurements are projection-noise limited then a 
Poissonian error model would be more accurate, but we 
shall see that this is nonetheless a good approximation. 

The main advantage of the latter formulation is that 
we can eliminate the explicit dependence on the linear 
coefficients a^^, hki, Cki and the noise variances ake by 
integration over suitable priors to obtain an explicit ex- 
pression for the probability of a particular model given 
the observed data d^^ that depends only on the six tran- 
sition frequencies ijJ, rather than the > 200 parameters in 
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the full model. Following standard Bayesian analysis 
we obtain 



P(a;|d)cx n 



1 



13(h| 



where the averages are defined by 



(13-Af)/2 



n=l 
1 " 



(16) 



(17a) 



(17b) 



The components hke-m are essentially the orthogonal pro- 
jections of the data onto a set of orthonormal basis vec- 
tors Hm{tn) 



N 



^ki:m 



Hm{tn)dkt-n- 



(18) 



The orthonormal basis vectors are derived from the (non- 
orthogonal) basis functions gm{i) defined above, evalu- 
ated at the respective sample times i„, via 



1 

■^mij^n) , / , ^m' rn9m,' (j^n) 



(19) 



where em'm is a 13 x 13 matrix whose columns e„i are 
the normalized eigenvectors — Gem — ame„i — of the 
13 X 13 matrix G = (Grn^ra^i) with 



N 

Grnim2 — ^ ^ ffmi (^w)ffm2 (^w) ■ 



(20) 



The objective is to find lj that maximizes P{u!\dke), or 
equivalently, the log-likelihood function 



13 -iV 



kJ=l 



(21) 

Note that N and (d^^) are constants, while hke indi- 
rectly depends on a; via the basis functions gm{t). It can 
be shown that the corresponding optimal coefficients are 
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FIG. 3: Estimated and actual values of the coefficients aki;m 
with estimated error-bars of System 1. When both A'' and 
A^e are large (top) the error-bars are nearly invisible, and 
the estimated and actual values are almost indistinguishable. 
When A'' and A^e are both small (bottom), the error-bars are 
significantly larger, mainly due to increased noise variances 
af.i\ yet the actual and estimated values for the coefficients 
are still almost indistinguishable. This suggests that the es- 
timated coefficients are in fact much more accurate than the 
uncertainty estimates suggest. 



We can similarly derive expressions for second moments 

13 

'\XM-rniXki-m2!~\Xkt,mij'\XM-m2! — ^kl 1 

m — 1 

(24) 

where ct^^ is the noise variance of the (fc, £)th signal, 
which can be approximated by its estimated expectation 



1 



>kil 



N -\h 



[N{dM) - 13(hw)] . (25) 



Note that for mi — m2 Eq. ( |24| is simply the variance 
of the parameter xi^g.^, which gives an estimate of the 
uncertainty Axke-m of the coefficient Xke-m 



A 



•^kt,m 



13 

E 



(26) 



afc^ = {{xM-i), (xke-s), • ■ ■ , (xki-ii)) , (22a) 

hki = {{Xki;2), {xM-i), . . . , {Xkt,12)) , (22b) 

Cki = {xki-iz), (22c) 

where {xkt-m) is shorthand notation for the expectation 
values E{xkl■,m\<^^ dki) of the linear coefficients of the ba- 
sis functions, given the optimal frequencies u) and the 
data dki- Furthermore ,22J , 



{Xkt-m) — 



m' — l 



am' 



(23) 



Fig. [3] shows that for a sufficiently large number of data 
points N and experiment repetitions per data point, N^, 
these uncertainties can be made very small indeed. For 
iV and/or small, the uncertainties are much larger, 
but simulations for our specific problem suggest that the 
estimated values are generally still very close to the ac- 
tual values even for small N and/or Ng, much closer than 
the uncertainty estimates would suggest. 



Although the log- likelihood function (21) depends ex- 
plicitly only on the six frequencies u) € rather than 
the full 214 model parameters, finding its (global) maxi- 
mum is not trivial as the log-likelihood is sharply peaked 
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FIG. 4: Power spectrum C{ijj) of System 1. Although the 
power spectrum is noisy, the log-plot of (7(0;) of the measured 
signals shown in Fig. [2] shows six well-defined peaks for LOm > 
in addition to the peak at a; = 0. The inset shows the 
filtered power spectrum C{uj) > Co, from which the six peaks 
ujm can easily be identified using standard peak detection. 
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FIG. 5: The power spectrum of System 12. The power Spec- 
trum C{u}) has only five peaks Um > in addition to the 
peak a.t Lo — 0. This could mean that the system has only five 
distinct transition frequencies, or that the measured signals 
are not sufficient to resolve two (closely spaced) transition 
frequencies. 



with many local extrema, and thus computationally ef- 
ficient gradient-based optimization algorithms are likely 
to get trapped in local extrema if the starting point is 
chosen randomly. An alternative is to use global search 
algorithms such as pattern search or genetic algorithms 
but these are computationally expensive and the results 
for our problem proved inaccurate. To circumvent this 
problem we adopt a combination strategy. 

We can first estimate the resonant frequencies by look- 
ing for peaks in the power spectra 



1 ^ 

ri = l 



(27) 



Using spectral filtering combined with a basic peak find- 
ing routine, we locate (up to) six peaks in the com- 
bined power spectrum 



(28) 



k,l=l 



as illustrated in Fig. |4j which are then used as input 
(i>(o) = (oj]^, . . . , cjg) to an optimization routine based 
on the BFGS quasi-Newton method with cubic line 
search [231 HH 113 US] to find the maximum of the log- 
likelihood (21 1. Although the discrete Fourier transform 



is not an optimal frequency estimator for a signal with 
multiple frequencies, it proved generally effective in pro- 
viding good starting values for the log-likelihood opti- 
mization routine, provided that the total sampling time 
(signal length) T was sufliciently long to resolve the res- 
onant peaks. 
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FIG. 6: (Color online) Log-likelihood on /i x Ji for Sys- 
tem 12 with five-peak power spectrum shown in Fig. [5] 
shows symmetry about y = x as \ogP{{uo\,u)2, ■ ■ = 
log P((a;2, tui, . . .)|dfcf) with twin peaks for y ^ x indicating 
that the most probable model on this subspace of the param- 
eter space is a six-frequency model. 



Since the frequency resolution of the power spectrum is 
limited by the signal length T, Acj — 1^, li there are two 
or more closely-spaced transition frequencies then it may 
not be possible to resolve six peaks in the power spec- 
trum without increasing the signal lengths significantly. 
But this is generally not necessary as we can improve 
the frequency resolution as follows. Suppose there are 
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Wl U)2 <^4 'i'S 1^6 


logP 




0.4293 0.8586 4.9983 5.4276 5.8569 
0.4291 0.8558 5.0046 5.4282 5.8604 


924.4486 
938.2960 




0.4235 0.4323 0.8558 5.0046 5.4282 5.8604 


943.3509 




0.4291 0.7631 0.8558 5.0046 5.4282 5.8604 
0.4291 0.8558 5.0046 5.1023 5.4282 5.8604 
0.4291 0.8558 5.0046 5.4282 5.5063 5.8604 
0.4291 0.8558 5.0046 5.4282 5.8604 5.9287 


938.3099 
938.2977 
938.2993 
938.2975 




0.4236 0.4322 0.8558 5.0046 5.4282 5.8604 





TABLE I: Log-likelihood for different five and six frequency 
models and actual transition frequencies for System 12 with 
five-peak power spectrum shown in Fig. |5] 



five identifiable peaks, Wi to oj^, in the power spectrum, 
as shown in the example in Fig. [5] Then we proceed as 
before, using the five peak frequencies in the power spec- 
trum as input u;^"^ for the optimization routine to find 
the most likely five- frequency model a;^*^. To ascertain 
whether there is a more probable six-frequency model 



we choose an interval Im about each 



,(*) 



TO = 1, 



and investigate the log- likelihood function (21) on the 
2D parameter space Im x /,„, keeping the other four fre- 
quencies fixed in each case. E.g., for to = 1 in the ex- 
ample above we find the maximum of \ogP{u>\Aki) for 
u) = (wi,ti;2,W2*\ti^3*\ a;4*'',u;5*'') with {uji,uj2) G and 

Ii — [io^i'^ — ^,uj^i^ + ^Y by calculating log P on a coarse 
2D grid, finding the maximum on the grid and using the 
resulting u; as a starting point for the BFGS optimization 
routine as before. A contour plot showing the maxima in 
the log-likelihood on Ji x Ii is shown in Fig. [6j 

We repeat this procedure for each to in turn. The re- 
sults, summarized in Table [Tj show that the six frequency 
model is most likely, more than the five-frequency 
model, and the other five six-frequency models. Indeed, 
the frequencies of the most likely six-frequency model 
are very close to the actual transition frequencies of the 
system simulated. However, the relative flatness of the 
peak corresponding to the global maximum of the log- 
likelihood function and the relatively small differences 
between the likelihood of the most likely model and the 
less likely models, suggests that more data would be de- 
sirable to improve the resolution of the parameter esti- 
mates, and our confidence that the model is indeed the 
correct choice. If there are fewer than five peaks in the 
power spectrum, the procedure described can be iterated 
to sequentially resolve peaks in the power spectrum until 
the most probable model has been found. 

To test the effectiveness and accuracy of this param- 
eter estimation technique, we test the method for 100 
randomly generated Hamiltonians, sampled at At ~ 0.1 
(arbitrary units) for different signal lengths T — (iV— 1) x 
M = 0.1 X 2'^ for d = 10, 11, 12, 13, 14 and different levels 
of projection noise, with the number of measurements per 
data point, £ {125, 250, 500, 1000}. The test Hamil- 
tonians have transition frequencies in the range of [0.3, 7], 



N\N, 


125 250 500 1000 


125 250 500 1000 


16,385 
8,193 
4,097 
2,049 
1,025 


0.093 0.094 0.094 0.094 
0.231 0.226 0.231 0.231 
0.432 0.432 0.432 0.432 
0.696 0.685 0.680 0.685 
1.646 1.650 1.646 1.650 


0.0002 0.0002 0.0001 0.0001 
0.0006 0.0006 0.0004 0.0003 
0.0018 0.0019 0.0009 0.0009 
0.0065 0.0040 0.0030 0.0024 
0.0272 0.0184 0.0085 0.0108 



TABLE IL The percentage relative errors {emax(i^'''')) (left) 
and {emax(<i'°''*)) (right) show that the log-likelihood opti- 
mization improves the accuracy of the frequency estimates by 
at least two orders of magnitude compared to the estimates 
obtained from the power spectrum. 





125 250 500 1000 


125 250 500 1000 


16,385 
8,193 
4,097 
2,049 
1,025 


0.068 0.068 0.068 0.068 
0.167 0.164 0.167 0.167 
0.327 0.327 0.327 0.327 
0.5100 0.493 0.493 0.493 
1.164 1.142 1.164 1.142 


0.0001 0.0001 0.0001 0.0001 
0.0005 0.0003 0.0002 0.0002 
0.0011 0.0012 0.0006 0.0005 
0.0035 0.0023 0.0019 0.0011 
0.0126 0.0089 0.0052 0.0036 



TABLE IIL The medians of percentage relative errors 
emax(<ij'') (left) and emax('i'°^*) (right) show the same accu- 
racy improvements of the log-likelihood estimates. Median 
errors lower than the averages indicate that the error distri- 
bution is peaked towards the origin. 



and include cases with very closely spaced transition fre- 
quencies, as shown in Fig. [7] To assess the quality of 
the models found, we calculate the transition frequen- 
cies Urn and corresponding parameters aki, '^m and Cki 
for each Hamiltonian, and consider the relative errors of 
the parameters identified from the noisy data with the 
parameter estimation technique described. 

Tables and |III| show the means and medians, respec- 
tively, over 100 systems, of the maximum relative error 
(in percent) 



,(a;(°)) = 100 X max 

rnGl....,6 



1 - 



(0) 



(29) 



of the estimated transition frequencies for each system, 
where are the exact transition frequencies. Com- 
parison of the errors for the initial frequency estimates 
obtained from the power spectrum, labeled ijJ^'^\ and 
the optimal values a;°P* obtained by maximizing the log- 
likelihood shows the optimized frequencies are generally 
about two orders of magnitude more accurate than the 
estimates obtained from the power spectrum. 

The linear coefficients akt,m, bke-m and Cke are then 
estimated from the maximization of Eq. (211 and 



from Eq.(22l. Taking the median of the relative errors 



emed (afc£;m) = 100% X median 



(30) 



where k, £ range from 1 to 4 and to = 1, . . . , 6, as a gen- 



eral measure of the quality of the fit. Table IV shows that 
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FIG. 7: (Color online) The transition frequency diagram for each of the 100 test systems shows that the transition frequencies 
range from 0.3 to 7, and there are six systems (12, 22, 34, 38, 73, 78) with two transition frequencies that differ by less than 
0.01 (circled), which are difficult to resolve, including one system (78) with two such cases. 





N\N, 


125 250 500 1000 


{emed(afcf;m)) 


16,385 
8,193 
4,097 
2,049 
1,025 


0.3825 0.2671 0.1912 0.1454 
0.5538 0.3598 0.2857 0.1923 
0.7711 0.5516 0.4075 0.2786 
1.0630 0.7940 0.5755 0.3762 
1.5817 1.1210 0.7880 0.5573 


(emod(bfe£;m)) 


16,385 
8,193 
4,097 
2,049 
1,025 


0.2417 0.1739 0.1174 0.0846 
0.3333 0.2519 0.1755 0.1144 
0.4860 0.3470 0.2394 0.1733 
0.6715 0.5098 0.3436 0.2485 
1.0194 0.7197 0.4691 0.3523 


{emcd(cfc^)) 


16,385 
8,193 
4,097 
2,049 
1,025 


0.0734 0.0525 0.0378 0.0279 
0.1002 0.0751 0.0538 0.0372 
0.1463 0.1037 0.0770 0.0518 
0.2007 0.1483 0.1148 0.0751 
0.2817 0.2258 0.1555 0.1047 




16,385 
8,193 
4,097 
2,049 
1,025 


0.0012 0.0006 0.0003 0.0001 
0.0971 0.0959 0.0953 0.0950 
0.2896 0.2873 0.2861 0.2855 
0.6763 0.6717 0.6692 0.6681 
1.4580 1.4487 1.4437 1.4414 



TABLE IV: Relative errors emcd(afcf;m), emod(?>fc<?;m), and 
emod(cfcf) (in %) and estimated error variances {(t^), averaged 
over 100 test systems for different signal length T — 0.1(A^ — 1) 
and number of experiment repetitions A'^e per data point. 



the average errors in the coefficients aki-m , bke-m and to a 
lesser extent Cke, are generally at least one order of mag- 
nitude larger than the error in the frequency estimates. 
Overall the quality is still good, however, with the (av- 
erage) errors ranging from a fraction of a percent to less 
than 2.5% for akt, and much less for Ckt, depending on 
the number of data points N and the accuracy of the data 
points determined by the number of experiment repeti- 
tions per data point, N^. Fig. |8] shows the distribution 
of the errors for both the least and greatest number of 
experiments. Apart from a few outliers, the distribution 
follows a roughly exponential form with most estimates 
being within a fraction of a percent of the true values, 
even for the least number of experimental samples. 

Table [n] shows that increasing Ng and thus the accu- 
racy of the data points does not improve the accuracy of 
the initial frequency estimates obtained from the power 
spectrum at all, while doubling N tends to reduce the er- 
ror by more than half. This is what we expect as once N^. 
is large enough to permit discrimination of the resonant 
peaks from the noise floor, little is gained by increas- 
ing Ne- Doubling does reduce the error for the opti- 
mized frequencies obtained from our Bayesian analysis, 
although if the accuracy of frequency estimates alone is 
considered, doubling the number of data points in prefer- 
able to doubling Ne- Increasing the accuracy (by dou- 
bling Ne) is more effective in reducing the errors in the 
coefficients a, b, c, but the contour plots in Fig. |9] show 
that the errors decrease faster with N, i.e., increasing the 
number of data points is generally still preferable. 





0.005 0.01 0.015 0.02 



Error % 

FIG. 8: Histogram of the relative % error for the test 100 
systems for samphng numbers of a) A'^ = 1025, A^e = 125 
b) A'^ = 16385, Ai'e = 1000. Inset graphs magnify the region 
around the origin showing the general distribution of errors 
which is roughly exponential. The numbers 66, 12, 78, 73 and 
77 in a) refer to outliers systems. 
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N 

FIG. 9: (Color online) Contour plots of the logj^Q mean (rela- 
tive) errors for the frequencies uj and coefficients a, b, c. The 
frequencies show the smallest errors (down to 10~®), whilst 
the a coefficients show errors up to a few percent (10"^) for 
the shortest signal lengths and greatest projection noise. 



tions {wia, u;24, '^14}, which must satisfy 



IV. HAMILTONIAN RECONSTRUCTION 



Wi3 = tJi2 + W23, (31a) 

UJ24: = (^23 + ^3i, (31b) 
UJu ^ UJ12 + LJ23 + ^SA- (31c) 



Once the frequencies u) and amphtudes a^.^, h^i and 
Cki have been extracted from the measured data using pa- 
rameter estimation, reconstructing the Hamiltonian (up 
to equivalence) requires at least two further steps: identi- 
fication of the resonant frequencies with transitions (/i, v) 
between eigenstates |^^) and l^i.) of the system, and com- 
putation of the parameters sm-i^ and Ake-^^i, in Eq. ( 10 ) 
from the coefficients a.ke, hki and Cke- For a four- level 
system we have three primary transitions {^12,^23,^34} 
between adjacent energy levels and three other transi- 



We identify the possible level structure (up to inver- 
sion) by examining the relationships between the frequen- 
cies. In the generic case, i.e., when there are six distinct 
transition frequencies, < uji < u;2 < • ■ • < 1^61 it fol- 



lows immediately from Eqs (31) that uiq — uju, and the 
primary transitions are {tt'i,W2,W6 — loi — uj2}- Closer 
inspection shows that there are 10 possible arrangements 
of the six transition frequencies as shown in Fig. |10| and 
the exact transition frequencies uj must satisfy AsU) = 
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FIG. 10; (Color online) Possible arrangements for a generic four-level system with six distinct transition frequencies. Not shown 
are the other five configurations which correspond to a refiection of the energies, which merely flips the above level structures. 



for one of the following matrices 



1110 
110-1 
110 

1110 
110-1 
10 10 






-1 




-1 



-1 




-1 






A. 



A,= 



1110 
10 1-1 
110 






-1 



110 10-1 
11-10 
10 1-10 



110 10-1 
11-10 
10 1-10 



(32a) 



Given the estimated frequencies 0;°^** the most likely case 
is that for which ||74sa;°P*||2 assumes its minimum, which 
should be close to 0, and significantly smaller than the 
errors for the other cases. A larger minimum error in- 
dicates and none of the possibilities is likely, suggesting 
that the system may not be a Hamiltonian four- level sys- 
tem. Similarly, if we have two cases for which the error 
the close to the minimum, this would be an indication 
that further data is required to resolve the ambiguity. 

Once the observed frequencies ujrn have been matched 
with actual transitions (/i, i^), we can associate the cor- 
responding coefficients ake,m, bki,m for m = 1, . . . , 6 with 
their respective transitions, i.e., we have aki-^i and bj^^.j^i, 
and determine the phase differences 



^ki-,iiv — arctan(6fc£.^^, afc£;^i,), 



(33) 



where arctan(6, a) is the four-quadrant arc tangent oib/a. 
If the estimated parameters are good, then the resulting 
Afcf;^^ should satisfy Akk;^iu ~ (mod 27r), Aki-f^u ~ 



-Ake-f^u (mod 27r), and 

A 
A 
A 



fe£;12 - 


\- Afc£;i3 


— Afcf;23 


= 


mod 


2tt 




1- Afc£;i4 


— Afc£;34 


= 


mod 




kt,12 - 


1- Afe£;14 


— Afc£;24 





mod 


2ti 



(34a) 
(34b) 
(34c) 



Due to the enforced symmetrization (15) of the coef- 
ficients a.ki and hki, the phase terms should satisfy 
Aki-^iv — —Aik-fiu- Minor violations of ([34| are to be 



expected, and can be mitigated, and the accuracy of the 
final reconstructed Hamiltonian improved by minimizing 
the constraint violations Hefe^Hl = X]s=i ^te-s^ where 

Cke-s = min{|a;fcf;s|, \xm-s - 27r|, \xke:s + 27r|}, (35) 
with Xfe£ = AA.ki and 



A 



110-10 
1 1 -1 
1 1 0-1 



(36) 



for k. 



1,- 



, 4 in a further refinement step, starting 
f^i, obtained from (33 1. This re- 



with the values for A 
finement tries to minimize the discrepancy between the 
estimated signal parameters and those expected from an 
underlying Hamiltonian model. It must be stressed, how- 
ever, that larger violations of the constraints are indica- 
tive of significant errors, which may even be exacerbated 
by such a refinement. In fact. Fig. 11 shows that there 
is a strong correlation between the maximum constraint 
violation prior to refinement 



max 

k4 



\ekf. 



(37) 



and the relative error of the final estimated Hamiltonian. 

Once the optimal values for Am;^^ have been found, 
we calculate the products 

Skht^Skt-v = akt.iiv cos(Afc£;^i,) + bkt-tiv sin(Afe£;^,^). (38) 
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FIG. 11: A scatterplot of the relative errors ||A/f||/||f/'|| of 
the estimated Hamihonian with llAi?]] as defined in (41 1 vs 



the maximum constraint violation (371 on a log-log scale for 
our 100 systems and 20 data sets per system (total of 2000 
data points) shows a strong correlation, suggesting that the 
maximum constraint violation (prior to refinement) is a good 
predictor of the accuracy of the estimated Hamiltonian. 



Labelling the RHS of the previous equation M^^.^,^ and 
defining the column vector s^^ and the 4x4 matrix Mm 



Ski = 



( SM-l\ 

Skt,3 



/Mki-ii ■■• Mke- 



M, 



ki 



ki-M 



\Mkt-Al 



Mkf: 



ki-Ai/ 



we can express Eqs ( 38 1 and (|9p) as follows 

SfcfSfc^ = Mm, slf,Si = Cki (39) 

for k, £ — 1, ... ,4:. To reconstruct the Hamiltonian ([lO| , 
we must determine the coefficients Ske;i, by solving ( |39[ ). 

Each Mm is a real symmetric matrix whose on- 
diagonal elements Mke-tii^, 1^ ^ ^, are determined by 
Eq. (38 1. The diagonal elements Mki-^fj,^ are unknown. 



However, we know that Mke should be a projector onto 
the ID space spanned by s^^, and the second equation in 
( 39 1 determines the norm of s^^ as well as the vector of 
diagonal elements (Mfcf.^^)^~Q. Thus, to determine the 
diagonal elements of Mke and the corresponding eigenvec- 
tor Sfe£, we note that a rank-1 projector H with matrix 
entries (gmn) must satisfy the condition 



dgn 



QmmQnn 



2 

9mn 



Vm, n. 



(40) 



Thus, given the off-diagonal elements of Mke, we choose 
the diagonal elements of Mm such as to minimize the 
norm of the error e — J2m n '^9mm and take s^^ to be 
the eigenvector corresponding to the eigenvalue of Mke 
closest to 1, normalized to ensure ||sfe£||2 = Cki- It is 
important to carefully choose the parameters for the op- 
timization here to ensure we find the diagonal elements 



N\N, 


125 


250 


500 


1000 




e[H) 1%5% 


e{H) 1%5% 


£{H) 1%5% 


£{H) 1%5% 


16,385 


0.40 11 1 


0.27 5 


0.18 2 


0.13 4 


8,193 


0.57 22 


0.41 8 


0.31 8 1 


0.19 4 


4,097 


0.87 41 5 


0.66 25 2 


0.41 15 1 


0.28 7 1 


2,049 


1.12 60 7 


0.91 45 6 


0.58 19 4 


0.44 12 2 


1,025 


1.81 81 13 


1.32 64 8 


0.84 34 5 


0.63 31 4 



TABLE V: Relative error £{H) = 100 x - H\\/\\H\\ 

of reconstructed Hamiltonian (with phase corrections) in %. 
Each table entry consists of three numbers: the median error 
(in %) and the number of systems (of 100) with relative error 
exceeding 1% and 5%, respectively. 



corresponding to the global minimum. Ideally, the resid- 
ual error e should be 10~^° or less. 

We implemented and tested the algorithm for our 100 
Hamiltonians. We were able to correctly identify the level 
structures for all but one case: system 73, which has two 
nearly identical transition frequencies with L02 = 1.8012 
and ujz = 1.8026, for N = 1025 data points sampled at 
Ne = 125, 250, and 500 experiment repetitions per data 
point. Even for this system, we were able to correctly 
identify the level structure by doubling the number of 
data points TV, with the exception of = 250 where at 
least N = 4097 data points were needed. Of course, in 
practice more data points would be required for such a 
system to be confident that the identification is correct, 
as explained earlier. 

To gauge the overall accuracy of the estimated Hamil- 
tonians we would like to compute the norm of the 
error ||AiJ|| — \\H'^^*- — iJ^'^'H, or the relative error 
||Ai7||/||iJ^'^'||, where we choose the operator norm here. 
However, calculating the norm of the error is complicated 
by the fact that we can only reconstruct the Hamiltonian 
up to the diagonal matrix D and energy inversion sym- 
metry. Thus we must compensate for the phases that are 
"unobservable" in our model by setting 

||Ai/|| = (41) 

with D = diag(l, 612, 61^, ^14), where 

5ie = phase(i/ff ) - phase(i?5=f ), £ = 2,3, 4, (42) 

and phase (i/J';'^') is the complex phase of the (1, Z) matrix 
element of i/^'^*, etc. Table |V shows the results of the 
percentage relative errors ||Ai/ for our 100 test 

systems, for different values of N and N^. Medians of 
the relative errors range from 0.13% for iV = 16385 and 
Ne = 1000 to 1.81% for TV = 1, 025 and = 125. 

V. CONTROL HAMILTONIAN TOMOGRAPHY 

We have seen that our procedure can characterize a 
single Hamiltonian up to a (physically irrelevant) global 
energy shift, and three relative phases (5i„ for n = 2, 3,4, 



12 



due to the freedom to redefine each of the measurement 
basis vectors by a U(l) phase minus an overall phase. If 
we can only measure the system in a fixed basis and pre- 
pare it in the measurement basis states, and the evolution 
is determined by a single fixed Hamiltonian, then wc have 
determined all observable parameters. However, for the 
system to be controllable, we require at least two (non- 
commuting) Hamiltonians, or more generally we must 
have the ability to modify the Hamiltonian by changing 
control parameters f, e.g., by applying external fields or 
varying applied gate voltages, etc. In this case we can 
still choose the phases S^^^J for one "reference" Hamilto- 
nian Hq = iJ(fo) as we wish, e.g., S^J = but the phases 
6^2 for all other Hamiltonians -ff(f) are now observable 
and thus relevant, and complete control Hamiltonian re- 
construction therefore requires that we identify them. 

To achieve this, note that if can initialize the system 
in the superposition state |<i>) — X]j=i measure 
the time-evolved state 



mt)) = Uf{m)=DlUr{t)Df\^) 



(43) 



with Uf{t) = exp[-itH{{)], Uf{t) = exp[-itH{i)] then 
Pi{t) = \{l\D\ij,{t)]Dfm'' = \{i\Uf{t)Dfm'' (44) 



shows that the phases (5[„ that determine Df = 
diag(l, e'*i2^ e'*i4) are now observable as Df no 
longer commutes with the initial state |<i>). As Uf{t) is 
fully determined by previous steps, if the initial state |<i>) 
is known, then the only unknown parameters in Eq. (44) 
are S^^ for n = 2,3,4. Given a set of measured values 
dik for pi(tk), we can determine the unknown parameters 
^in by minimizing the least-squares error 



e = 



El 



(45) 



where Pi = {pi{ta), . . . ,pi{tK)) and = {dj>o, dix) 
for 1= 1,2,3,4. An explicit expression for pi{t) derived 
in Appendix \K\ shows that we can in principle determine 
all the phases if the initial state satisfies aj ^ for all 
j [29 . Moreover, it is advantageous to choose a balanced 
initial state, |ajp ~ \ for all j, if possible, to maximize 
signal to noise ratios. 

To prepare such an initial state, we can use the ref- 
erence Hamiltonian Hq. Unless the reference Hamilto- 
nian is such that one or more of the measurement ba- 
sis states are completely decoupled from state |1), it 
is almost certain that the time-evolved state \^\(t)) = 
t^o(t)|l) = Y.]=iaj{t)\3) with Uo{t) = exp{-itHo) will 
satisfy aj{t) for all j for at least some t > 0. Thus, 
having characterized the Hamiltonians H = H{i) for 
different control settings f up to the phases (J^^, all we 
need to do is to select a suitable reference Hamiltonian 
Hq = _ff (fg), and find a time such that the time-evolved 




FIG. 12: (Color online) Evolution of populations aj{t)\ un- 
der reference Hamiltonian (system 5) and error J^^. \aj\^ — | 
from ideal balanced initial state. We selected the second mini- 
mum (which is the global minimum for < t < 10) a.tt = 5.34 
as initial evolution time for the 5 estimation step. 
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state |<I>i(i*)) satisfies \aj{t^ 



This is generally not 



FIG. 13: (Color online) Median relative error of estimated 
Hamiltonian after 5 estimation as a function of signal lengths 
A'' and projection noise level A^e of the measured traces pe{tk)- 
In all cases the most accurate estimates for H from the pre- 
vious step, i.e. A'^ — 16385 &i Ne — 1000, were used. 



difficult. For instance, we randomly choose the Hamil- 
tonian for test system 5 as our reference Hamiltonian. 
Fig. [12] shows that there are several times t e [0, 10] at 
which the populations |aj(t)p of all levels (in the mea- 
surement basis) are approximately equal. We pick one 
of these times — 5.34, set |<i>) — C/o(io)|l)i obtain 
the measurement traces pi{tk) as follows: 

1. Initialize system in measurement basis state |1). 

2. Let it evolve under Hamiltonian Hq for time i*. 

3. Change control settings to f and let system evolve 
for t time units under Hamiltonian Hf. 

4. Perform measurement outcome £ — 1,2,3,4. 
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N\N, 


125 
SiH) 1%5% 


250 

£{H) 1%5% 


500 

8iH) 1%5% 


1000 

8{H) 1%5% 


16,385 


0.89 40 2 
2.12 78 12 


0.70 27 1 
1.54 75 9 


0.56 14 1 
1.00 50 6 


0.45 6 1 
0.76 35 1 


8,193 


1.09 54 5 
2.60 89 24 


0.84 35 
2.21 86 14 


0.69 28 1 
1.66 77 10 


0.53 17 1 
1.08 52 5 


4,097 


1.48 68 7 
3.47 96 41 


1.12 58 5 

3.13 94 27 


0.91 43 5 
2.20 87 13 


0.61 26 2 
1.45 65 7 


2,049 


2.24 88 15 
6.06 98 67 


1.45 74 7 
3.91 93 37 


1.12 55 7 
2.85 89 25 


0.78 37 4 
2.25 78 16 


1,025 


3.14 95 29 
8.36 98 76 


2.44 90 18 
5.92 96 59 


1.60 80 8 
4.50 95 48 


1.22 59 6 
3.00 88 32 



TABLE VI: Relative error £{H) = 100 x - H\\/\\H\\ 

of reconstructed Hamiltonian with estimated phases Sin in % 
(no phase corrections). As before, each table entry consists 
of three numbers: the median error (in %) and the number 
of systems (of 100) with relative error exceeding 1% and 5%, 
respectively. The first row in each box are the estimates ob- 
tained for signals pi{tk) of length A'' — 50, the second row 
the estimates obtained for signals of length N = 200, in both 
cases sampled at A^e = 5000. 



As before we repeat this experiment N^. times for a fixed 
t to estimate pe{tk) (number of times the outcome was £ 
divided by Nf,), and then repeat for different times tk to 
obtain estimates for pe(tk). 

We tested the phase estimation procedure for the 
estimated Hamiltonians obtained in the previous step. 
For each of the 100 systems we first generated (simu- 
lated) measurement signals for pi(tk) of varying length 
T = {N — l)At and levels of projection noise N^. The 
number of points ranged from A'^ — 1 = 25 to 1000 data 
points, sampled at At = 0.1 fixed as before; the mea- 
surement repetitions from 1000 to 5000. In the re- 
construction of the phases, we only assume we know the 
estimated Hq, hence the estimated |<i>(0)), and the es- 
timated Hf, determined in Section |TV] While the most 
accurate estimates for the frequency and linear coefficient 
estimation step (step 1) were obtained for the longest sig- 
nals (A^ — 16, 385), we find that the accuracy of the phase 
estimation step peaks at around A^ sa 50, and that longer 



signals are in fact highly detrimental (Fig. 13 1. This may 
seem very surprising at first but can be at least partly 
explained by the fact that even small inaccuracies in the 
initial estimates, especially for the frequencies, will accu- 
mulate over time and increase the discrepancy between 
the projected evolution of the system based on our Hamil- 
tonian estimates and the true evolution. 

Based on these results we settled for signals of length 
A^ — 1 = 50 with A"e = 5000 measurement repetitions 
per data point for the final phase estimation step. For 
each of the 2000 estimated Hamiltonians H obtained in 
the first step — corresponding to the 100 different test 
systems, as well as four levels of projection noise A'e G 
{125, 250, 500, 1000} and five signal length T = {N-l)At 



for N-1 e {2^0, 2", 2^2, 2^3, 2"} with At = 0.1 fixed, 
each — we estimated the phases (5i„, and used the re- 
sults to reconstruct the total Hamiltonian H = HD. 
Table IVTl shows the results in terms of the median of rel- 
ative errors. For comparison we include the Table the 
results obtained had signals of length A^ = 200 been 
used instead. Comparison of the numbers clearly shows 
that longer signals are detrimental for the phase esti- 
mation step. In addition to substantially decreased ac- 
curacy, longer signals also slowed down the numerical 
optimization, making it more difficult for the routine to 
find the global minimum. In view of the complicated de- 
pendence of p£{tk) (see appendix [A| on the parameters 
Sin, n — 2,3,4, we initially explored population-based 
(global) optimization strategies, especially evolutionary 
algorithms, but found that it was substantially slower 
and far less effective in finding the global minimum of 
the error Eq. (45 1 than a gradient-based (BFGS-type) lo- 
cal optimization algorithm. In fact, for short signals the 
local optimization routine generally succeeded in finding 
the global minimum in a single run, starting with a ran- 
dom guess for d = (<^i2: ^13: ^14)1 although the optimiza- 
tion was repeated with several different initial guesses to 
increase the probability that we had indeed found the 
(globally) best value for 5. 



VI. CONCLUDING DISCUSSION 

We have presented a method for characterizing the 
Hamiltonian and its dependence on external control pa- 
rameters, which is a pre-requisite for Hamiltonian Engi- 
neering and coherent control of the system's evolution, 
for a generic two-qubit system, assuming only the ability 
of preparation and measurement in a fixed basis. Anal- 
ysis of simulated measurement data shows that the task 
of estimating the parameters from the complex, noisy 
measurement signals with multiple frequencies, and re- 
constructing the Hamiltonian is very challenging, and 
requires a carefully designed multi-step approach, com- 
bining spectral analysis, Bayesian analysis and several 
carefully designed optimization steps to reconstruct the 
energy level structure and matrix representation of the 
Hamiltonian. In the absence of any control, the Hamilto- 
nian can only be reconstructed up to three phases, due to 
the freedom to redefine the measurement basis by U(l) 
phase rotations. This symmetry can be broken if the 
system can be prepared initially in a suitable superposi- 
tion state, and we exploit this fact to achive full control 
Hamiltonian tomography in a simple two step procedure. 

The Bayesian analysis assumes a Gaussian noise pro- 
file which, though not strictly accurate, works well, es- 
pecially in the large A^e limit. Any significant deviations 
from Gaussian noise (e.g. Poissonian statistics for small 
Ne for p 0,1) will tend to make the log-likelihood 
estimates worse, and thus our estimates of the confi- 
dence that the model fits the data are conservative P^ . 
More accurate error estimates could be obtained using 
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Bayesian analysis with a Poissonian noise model, though 
our results show that even a Gaussian noise model re- 
sults in a huge improvement of two orders of magnitude 
or more in the accuracy of the frequency estimates, com- 
pared to estimates obtained from simple spectral analy- 
sis. This turned out to be crucial for successful Hamil- 
tonian reconstruction. The frequency estimates obtained 
from the power spectrum combined with a simple least- 
squares error minimization to find the optimal spectral 
amplitudes proved to be too inaccurate for Hamiltonian 
reconstruction, leading to inconsistent equation systems 
and significant errors, and any attempt to obtain esti- 
mates of the parameters by direct minimization of the 
least-squares error of the measurement signals and the 
expected signals resulted in reconstructed Hamiltonians 
that were little better than random for our test systems. 

Though we have implicitly assumed a Hamiltonian 
model, i.e., that incoherent effects will be negligible on 
the time scales of interest, any significant deviation from 
the assumed model, e.g., significant decoherence or cou- 
pling to additional states outside the two-qubit subspace 
would result in low likelihoods of the chosen (four-level) 
Hamiltonian model. Such effects can easily be incorpo- 
rated into the analysis by changing the basis functions, 
e.g., using damped exponentials instead of sinusoids or 
including additional states, which we will consider in fur- 
ther work. Furthermore, any prior information about the 
structure of the Hamiltonian can be incorporated to make 
the Bayesian analysis more efficient. Thus, the method 
lends itself to adaptive protocols, as we can adaptively 
sample the system until certain targets for the likelihood 
or error estimates are met, ensuring that we perform 
enough measurements to get accurate estimates but no 
more than necessary. |30j This is especially important as 
the number of measurements required will vary depend- 
ing on the system. For instance, for a system with well 
spaced transition frequencies, a sharply peaked likelihood 
function with a clearly identifiable global maximum can 
be obtained with much less data than for a system with 
two almost degenerate transition frequencies. 



For control Hamiltonian tomography, the small but 
non-zero inaccuracies in the initial estimation step lead 
to an optimum sampling time for the second step due to 
divergence of the model from the true system behavior 
at longer times. In principle, it should be possible to use 
this divergence to improve the initial estimates of the 
Hamiltonians, and exploring such refinements could be 
an interesting avenue for future research. Errors in the 
second step decreased with increased signal to noise ratio 
(increasing N^), as the estimate of the phase parameters 
does not depend on the signal length, unlike frequency 
resolution. It would also be interesting to investigate 
the accumulation of errors in this multi-step estimation, 
especially how uncertainties in prior steps affect the ac- 
curacy of the Bayesian estimation in subsequent stages. 
Finally, in this paper we have dealt with the generic case. 
When the Hamiltonian has exact degeneracies then the 
measurement signals will contain fewer than six frequen- 
cies. In this case, the level structure reconstruction be- 
comes harder as the number of special sub-cases increases 
and we may not be able to uniquely identify the Hamilto- 
nian. Although the set of Hamiltonians with exact degen- 
eracies is of measure zero, further study of these special 
cases may be of interest as one may want to specifically 
engineer Hamiltonians with such level structures. 
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APPENDIX A: MEASURED PROBABILITIES 

If the system is initialized in the generic superposition 
state |<i>) = ^^"i measured after evolving for 

t time units under the Hamiltonian Hb = HbD, then 
the general expression for the probability pi(t) of mea- 
surement outcome £ is 
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